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Abstract. Dense packings composed of irregular polyhedral particles are investigated by numerical simulations under 
quasistatic triaxial compression. The Contact Dynamics method is used for this investigation with 40 000 particles. The effect 
of particle shape is analyzed by comparing this packing with a packing of similar particle size distribution but with spherical 
particles. We analyze the origin of the higher shear strength of the polyhedra packing by considering various anisotropy 
parameters characterizing the microstructure and force transmission. Remarkably, we find that the polyhedra packing has 
a lower fabric anisotropy in terms of branch vectors (joining the particle centers) than the sphere packing. In contrast, 
the polyhedra packing shows a much higher force anisotropy which is at the origin of its higher shear strength. The force 
anisotropy in the polyhedra packing is shown to be related to the formation of face-face contacts. In particular, most face-face 
contacts belong to strong force chains along the major principal stress direction whereas vertex-face and edge-edge contacts 
are correlated with weak forces and oriented on average along the minor principal stress direction in steady shearing. 
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INTRODUCTION 

During the last two decades, granular media composed of 
circular particles (in 2D) and spheres (in 3D) have been 
a subject of systematic research. In particular, various 
microscopic features such as fabric anisotropy |1], force 
transmission f^l, (sl, 13, ^sl, ^6] and friction mobilization 
JtI] have been analyzed numerically and experimentally. 
Hence, an emerging issue today is how robust these 
findings are with respect toparticle properties such as 
shape and size distribution 18 91 Mon. 

In this paper, we study numerically granular materials 
composed of polyhedral particles. We use the contact 
dynamics method to simulate the slow shear of these 
materials in comparison to sphere packings with similar 
particle size distribution. The facetted shapes give rise 
to a rich microstructure where the particles touch at 
their faces, edges and vertices. We analyze the fabric 
and force anisotropics and their link with the stress- 
strain behavior. We show that face-face contacts play a 
major role in force transmission and statics of polyhedra 
by accommodating long force chains that are basically 
unstable in a packing composed of spheres. 



NUMERICAL PROCEDURES 

The simulations were carried out by means of the contact 
dynamics (CD) method with irregular polyhedra parti- 
cles 111, 12]. We used LMGC90 which is a multipurpose 
software developed in our laboratory, capable of model- 




FIGURE 1. Examples of polyhedra used in the simulations. 



ing a collection of deformable or undeformable particles 
of various shapes by different algorithms lisll . 

We generate two numerical samples. The first sam- 
ple (SI) is composed of 36933 polyhedra. The particle 
shapes are taken from a library of 1000 digitalized ballast 
grains provided by the French Railway Company SNCF. 
Fig. [T] shows several examples of the polyhedral particles 
used in the simulations. We used the following size dis- 
tribution: 50% of diameter dmin = 2.5 cm, 34% of diam- 
eter 3.75 cm, 16% of diameter dmax = 5 cm, where dmin 
is defined as two times the largest distance between the 
bary center and the vertices of the particle. The second 
sample (S2) is composed of 19998 spheres with exactly 
the same size distribution as in SI. Fig. [2l shows a snap- 
shot of the sample SI in equilibrium state after deposi- 
tion and isotropic compression under a constant stress of 
Oo = 10"^ Pa in a rectangular box at zero gravity. 

The coefficient of friction is 0.5 between the particles 
and with the walls. The initial value of the solid fraction 
is p :^ 0.6 in both samples. Both samples have a nearly 




FIGURE 2. A snapshot of the packing SI (polyhedra). The 
walls are not shown. 



square bottom of side such that L ^ / and an aspect ratio 
where H is the height. The isotropic samples 
were subjected to vertical compression by downward 
displacement of the top wall at a constant velocity of 
1 cm/s for a constant confining stress (To acting on the 
lateral walls. 



SHEAR STRENGTH 

In this section, we compare the stress-strain behavior 
between the packings of polyhedra (packing SI) and 
spheres (packing S2). For the estimation of the stress 
tensor, we use the internal moment tensor of each 
particle i defined by lll4ll : 

where is the a component of the force exerted on 
particle i at the contact c, is the j3 component of the 
position vector of the same contact c, and the summation 
runs over all contact neighbors of particle i. 

It can be shown that the internal moment of a collec- 
tion of rigid particles is the sum of the internal moments 
of individual particles [14]. The stress tensor G for a 
packing of volume V is simply given by 

where is the branch vector joining the centers of the 
two touching particles at the contact c. 

Under triaxial conditions with vertical compression, 
we have a\ > (72 = (73, where the Ga are the stress 
principal values. We extract the mean stress p = (cJi + 
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FIGURE 3. The normalized shear stress ^//? as a function of 
shear strain £q for the polyhedra packing S 1 and sphere packing 
S2. 
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FIGURE 4. Geometry of a contact between two polyhedra. 

<72 + 0"3)/3, and the stress deviator ^ = (cTi — (73)/3. For 
our system of perfectly rigid particles, the stress state is 
characterized by the mean stress p and the normalized 
shear stress q/p. 

The cumulative strain components Ea are defined by 

where Hq, Iq and Lq are the initial height, width and 
length of the simulation box, respectively, and AH = 
Hq — H, Al = Iq — I and AL = Lq — L are the cor- 
responding cumulative displacements. The cumulative 
shear strain is defined by £q = £1 — £0,. 

Figure [3] displays the evolution of q/p for the pack- 
ings SI and S2 as a function of £q. For both packings, we 
observe a classical behavior characterized by a harden- 
ing behavior followed by (slight) softening and a stress 
plateau. The higher level of q/p for the polyhedra pack- 
ing reflects the organization of the micro structure and the 
features of force transmission that we analyze below. 

GEOMETRICAL ANISOTROPY 

For the analyses that will be discussed below, we intro- 
duce the local frame {n' ^t') where n' is the unit vector 
along the branch / and t' is an orthonormal unit vector ; 
figure m We set 

l = ln', (4) 

where I is the length of the branch vector. 

We define the angular averages associated with the 
branch vectors Let J2/(^1) be the set of branch vectors 
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FIGURE 5. Polar representation of the probability density 
function Pq^ and (^(^1)) for the samples SI and S2 in the 
residual state. 



pointing in the direction ^1 = (0, 0) up to a solid angle 
dQ. and A/c(^^) its cardinal. The angular averages are 
defined as follows: 



1 



(5) 



where Nc is the total number of contacts, and 1^ the 
actual values of branch vector length, for contact c, re- 
spectively. 

Under the axisymmetric conditions of our simulations, 
these two functions are independent of 0. Fig. [5] displays 
a polar representation of these functions in the -plane 
for polyhedra (SI) and spheres (S2) at Eq = 0.3. We 
observe an anisotropic behavior of the unit inter-center 
vector Pq{6) in both cases. A weak anisotropy of branch 
vector can be seen for S 1 . The magnitude of anisotropy 
is larger for spheres compared to polyhedra except for 
{i){9) which is weakly anisotropic for polyhedra. 

The simple shapes of the above functions can be 
approximated by harmonic approximation. Considering 
only the functions compatible with the symmetries of the 
problem (independent with respect to (j) and TT-periodic 
as a function of 0) we have 



Pa{0) = -^-{l+a[3cos^e-l]}, 
471 



4{ 1+^/ [3cos^e-l] } 



(6) 
(7) 



where a, ai are the anisotropy parameters, £o is the mean 
branch vector length. The probability density function 
/h(0) is normalized to 1. The harmonic fits are shown 
in figure [5] for the two functions in the critical state. 

The evolution of the anisotropics with Eq are displayed 
in Fig. [6] for SI and S2. We see that a is systematically 
larger for spheres than for polyhedra. The branch vec- 
tor length anisotropy ai is negligible for spheres. The 
low anisotropy of the polyhedra packing results from a 
particular organization of the force network in correla- 
tion with the orientations of each contacts (edge-to-face, 
vertex-to-edge...) in the packing 




FIGURE 6. Evolution of anisotropics a, ai with Sq for pack- 
ings SI and S2. 



FORCE ANISOTROPY 

We consider the components of the contact force by 

f = fn'n'^ft't\ (8) 

We refer to and as radial and orthoradial com- 
ponents of the contact force. As the angular orientation 
of the branch vector, we distinguish the angular distri- 
butions of radial forces (/„/)(^l) and orthoradial forces 
{fft ) {Q). These two functions can be expanded on a base 
of spherical harmonics. At leading order, we have IslfTsIl 



{ft'){0) 



fo{ 1 [3cos^ 
fo aff sin 20, 



1]}, 



(9) 



where a^/ and af/ are the anisotropy parameters, and fo 
the mean force. The anisotropics a^f and a^f are plotted 
in figure[7]as a function of £q. The radial force anisotropy 
a^f increases as the fabric anisotropy, and tends to a 
plateau. But, in contrast to fabric anisotropy, its value 
is higher for polyhedra than for spheres. This means 
that the large force anisotropy is correlated with particle 
shape rather than with fabric anisotropy. The orthoradial 
force anisotropy a^f has a similar behavior except that it 
takes considerably higher values in the case of polyhedra 
compared to spheres due to large friction developed by 
face to face contacts 

The anisotropics a, ai, a^f and a^f are interesting de- 
scriptors of granular micro structure and force transmis- 
sion as they underlie the shear stress. Indeed, it can be 
shown that the general expression of the stress tensor Eq. 




FIGURE 7. Evolution of anisotropies af with Eq for 
packings SI and S2. 
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FIGURE 8. The normalized shear stress q/p as a function 
of shear strain £q for the packings SI and S2 both from direct 
simulation data and theoretical prediction of Eq. ( [TOl ). 



^ under some approximations leads to the following 
simple "stress-force-fabric" relation iHlsl 



P 



(10) 



As we see in Fig. [8l our simulation data are in quan- 
titative agreement with this relation both for spheres 
and polyhedra, all along the shear. A remarkable conse- 
quence of Eq. ([TOl) is to reveal that the fabric anisotropy 
provides a major contribution to shear stress in the sphere 
packing whereas the force anisotropies are more impor- 
tant for shear stress in the polyhedra packing. 



CONCLUSION 

The objective of this paper was to isolate the effect of 
particle shape with respect to shear strength in 3D gran- 
ular media by comparing two similar packings with dif- 
ferent particle shapes. A novel finding of this work is 
that the origin of enhanced shear strength in a poly- 
hedra packing compared to a sphere packing lies in 
force anisotropy induced by particle shape. The fabric 
anisotropy associated with the network of branch vec- 
tors is lower in the polyhedra packing. In other words, 
the force anisotropy, partially underlying shear strength, 
is mainly controlled by the fabric anisotropy in a sphere 
packing. This mechanism breaks down to some extent 
in a packing of polyhedra where force anisotropy results 
mainly from the "facetted" particle shape. 
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